arXiv:1503.03894vl [cond-mat.mtrl-sci] 12 Mar 2015 


Atomistic simulations of amorphous polymers 
in the cloud with PolymerModeler 


Benjamin P. Haley\ Chunyu Li^, Nathaniel Wilson^, Eugenio 
Jaramillo^, and Alejandro Strachan*^ 

^Research Computing 

^School of Materials Engineering and Birck Nanotechnology 
Center, Purdue University, West Lafayette, Indiana 47907 
^Department of Biology and Chemistry, Texas A&M International 
University, Laredo, TX 78041 

March 16, 2015 


Abstract 

Molecular dynamics (MD) simulations enable the description of ma¬ 
terial properties and processes with atomistic detail by numerically solv¬ 
ing the time evolution of every atom in the system. We introduce Poly¬ 
merModeler, a general-purpose, online simulation tool to build atomistic 
structures of amorphous polymers and perform MD simulations on the re¬ 
sulting configurations to predict their thermo-mechanical properties. Poly¬ 
merModeler is available, free of charge, at nanoHUB.org via an interactive 
web interface, and the actual simulations are performed in the cloud us¬ 
ing nanoHUB.org resources. Starting from the specification of one or more 
monomers PolymerModeler builds the polymer chains into a simulation cell 
with periodic boundary conditions at the desired density. Monomers are 
added sequentially using the continuous configuration bias direct Monte 
Carlo method, and copolymers can be created describing the desired se¬ 
quence of monomers. PolymerModeler also enables users to perform MD 
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simulations on the structures created by the builder using the Dreiding force 
field and Gasteiger partial atomic charges. We describe the force field im- 
plemenfafion and fhe various opfions for fhe MD simulations fhaf use fhe 
LAMMPS simulafor. PolymerModeler MD simulations for a PMMA sam¬ 
ple show good sfrucfural agreemenf wifh experimenfs and are in good agree- 
menf wifh simulations obfained wifh commercial soflware. 


1 Introduction 

Atomistic simulations are a powerful tool to predict materials properties and un- 
eover the fundamental meehanisms that govern them from first prineiples; these 
teehniques are expeeted to play an inereasingly important role in the design of new 
materials with tailored properties. |[Ill2l Among many other eontributions, molee- 
ular dynamies (MD) simulations have provided invaluable insight into the re¬ 
sponse of materials under extreme eonditions|l3l|4l|, size effeets on plastie deform- 
ation[|5l. and polymer dynamies MD simulations are ubiquitous in seienee 
and engineering due, in part, to the eontinuing inerease in eomputing power and 
the availability of efficient, parallel, and scalable eodes that ean make effective 
use of current supereomputers [r7l[^l9l [T0l[TT]| . 

Even though some of these simulation eodes, sueh as ESPResSo lfT^ . 

GROM ACS IfT^ . and AMBER IfT4l are designed to simulate soft matter systems, 
the simulation of amorphous polymers remains ehallenging due to the diffieulties 
involved in generating aeeurate moleeular struetures to be used as initial eondi- 
tions in MD simulations. Unlike erystalline materials where initial atomie posi¬ 
tions ean be easily generated by replieating a unit eell, or in atomic amorphous 
solids where equilibrated struetures ean be obtained from the melt [fT5l [T^ . the 
generation of amorphous thermoset or thermoplastie polymer struetures requires 
speeial-purpose algorithms [fT7l[T^[T9ll20ll^l22ll2^ . Chain dynamies in linear- 
ehain polymers involve timeseales well beyond those aehievable with MD simu¬ 
lations, espeeially for large molecular-weight cases lf24ll . and if the initial strueture 
does not exhibit eorreet statistieal properties the simulations will produee inaeeu- 
rate results. 

Some eommereial software paekages, sueh as MAPS by Seienomics and Ma¬ 
terials Studio by Aeeelrys, Ine, offer the eapability to build atomistie struetures of 
eondensed-phase amorphous polymers using a simple GUI, paek the ehains into a 
simulation volume, and perform MD simulations. Open souree eodes designed to 
build amorphous polymers inelude Polymatie ||22ll2^, MCCCS Towhee lf25ll2^ . 
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PACKMQL EtI . AMBER IfT4l . and GRQMACS lfT3ll . These eodes requires users 
to ereate input files in domain-speeifie mini-languages, in addition to eompiling 
and running the program in order to paek ehains in preparation for MD simula¬ 
tions. The ECCE [|2^ tool offers a GET for building systems that may be simu¬ 
lated using ab initio paekages sueh as Gaussian [[29l but not with general purpose 
MD eodes. Thus, we designed and developed PolymerModeler to provide a easy- 
to-use, general-purpose, freely aeeessible tool to build amorphous thermoplastie 
struetures and perform MD simulations. We expeet the tool to be useful both for 
edueation and in researeh: setting simulations up is simple enough to be used in 
the elassroom, all that is required is an internet eonneetion, yet, the eode is pow¬ 
erful enough for researeh and ean also be used in eonjunetion with state-of-the-art 
teehniques to relax polymer struetures and prediet their properties. 

The PolymerModeler tool has been deployed on nanoHETB.org and is freely 
available for online simulations via the URL http://nanohub.org/resourees/polymod 
The only requirement is a nanoHUB aeeount that ean be obtained after a free reg¬ 
istration. nanoHUB [f30l makes simulation tools aeeessible through a standard 
web browser or an iPad app, removing the need for users to download, eonfigure, 
and install software and learn seripts or domain-speeifie languages. The state of 
a simulation is preserved even after elosing the browser in whieh it is running or 
logging out of nanoHUB so that the user ean return to the simulation later. Users 
ean provide feedbaek and ask questions in nanoHUB.org. This information has 
guided the development of PolymerModeler, and we will eontinue listening to the 
desires of the user eommunity for future development. 

The PolymerModeler tool features a builder of amorphous, linear-ehain ther- 
moplastie polymers, detailed in Seetion[2l As deseribed in Seetion[3l the atomistie 
system eonstrueted by the builder ean be used as initial eonditions for MD simu¬ 
lations using LAMMPS OTl to eompute thermo-meehanieal properties of the 
polymer. All simulations run on nanoHUB eomputing resourees ineluding the 
possibility of running the MD in parallel over multiple proeessors. Conelusions 
are drawn in Seetion[5l 


2 Amorphous polymer builder 

The proeess of building a eondensed-phase, amorphous, linear-ehain polymer in 
PolymerModeler ean be divided in three main steps: i) the speeifieation of one 
or multiple monomers and their statistieal arrangement into ehains, ineluding the 
ehain length and number of ehains to paek in the simulation eell, ii) the deter- 
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mination of the simulation cell into which the chains will be placed; and iii) the 
specification of rules for the sequential placement of monomers into chains of 
the desired length with conformation given by a set of dihedral angles between 
monomers. In PolymerModeler this last step is performed using the continuous 
configurational bias Monte Carlo algorithm IfTTl . The following subsections de¬ 
scribe these steps including details of the algorithms, approximations, and their 
implementation after which results of the builder are presented. Once the chains 
are built, PolymerModeler performs a series of relaxations to remove bad con¬ 
tacts and reduce entanglements after which molecular dynamics can be performed 
to predict thermo-mechanical properties; these relaxation and MD steps are de¬ 
scribed in Section [3l 

2.1 Monomers and their arrangement 

The first step in building amorphous thermoplastic polymers is specifying the re¬ 
peating unit which will be used to construct chains, and users can choose between 
a set of pre-built monomers or create/upload their own. The PolymerModeler tool 
accepts monomer specifications in the commonly used PDB (protein data bank) 
and XYZ file formats. Files of these types are available from a wide variety of 
online databases and may be uploaded directly to the PolymerModeler tool with 
no modification. The user must also specify which atom in the monomer is the 
head atom and which is the tail. The builder constructs a chain by adding one 
monomer unit at a time. Once a new monomer is added, the tail atom of the 
previous monomer on the chain is removed and so is the head atom of the new 
monomer. The final step is the creation of a bond between the backbone atoms 
originally attached to the head and tail atoms just removed. The head and tail 
atoms must each be connected to the molecule by a single bond. 

Internally the PolymerModeler tool converts a PDB or XYZ monomer to a 
z-matrix (internal coordinates format) representation that details the bonds, bond 
lengths, bond angles, and torsion (dihedral) angles between atoms. In this repre¬ 
sentation, the head atom is listed first, followed by all the backbone atoms to the 
tail, then all other atoms in the monomer. The PolymerModeler tool also accepts 
monomer specifications in this z-matrix format. The box describes the z-matrix 
format and shows an example for a PMMA monomer. 
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Figure 1: Internal coordinate representation of a PMMA repeated unit 


The z-matrix format 

The first line specifies the number of backbone atoms in the monomer, 
which are the leading entries in the z-matrix. The first column of the 
z-matrix gives information about atomic species in terms of Dreiding 
atomic types (discussed in the following paragraph). From the second 
atom on, the second column of the z-matrix indicates the index of a 
previously defined atom to which the current atom (the atom specified 
by the current line) is bonded; the third column is then the bond length, 
in Angstroms, between these two atoms. For the third atom on, an angle 
needs to be specified to fully determine its internal position, so two 
additional columns are needed indicating a third atom (column four) 
and the corresponding bond angle, in degrees, is in the fifth column. 

For the remaining atoms a dihedral angle needs to be specified; the sixth 
column gives the index to a fourth atom required to specify a dihedral 
angle (defined between the plane determined by atoms on columns 6, 4 
and 2 with that determined by atoms on columns 1, 2 and 4). The torsion 
angle between these planes is given, in degrees, in the final column. 
Pre-built z-matrix representations of several monomers are available in 
PolymerModeler. 

Once the list of monomers are specified, the PolymerModeler builder accepts 
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a list of monomer arrangements, eonsisting of repeating patterns or probabili¬ 
ties to ereate random struetures. A pattern arrangement speeifies an order in 
whieh monomers are added to a ehain, repeating until the desired ehain length 
is reaehed. A probability arrangement assigns relative weights to the defined 
monomers, whieh are stoehastieally added to a ehain aeeording to these weights. 
Multiple arrangements (patterns or probabilities) ean be defined, with relative 
weights, to support the growth of multiple ehain eonfigurations in the same sys¬ 
tem. The builder ehooses a pattern for eaeh ehain added to the system, aeeording 
the the weights assigned to the patterns, and eonstruets the entire ehain aeeording 
to the seleeted pattern. In this manner the builder allows the full speeifieation of 
stereoehemistry. The supplementary material provides examples of how these fea¬ 
tures ean be used to ereate isotaetie, syndiotaetie, and ataetie PMMA struetures. 

2.2 Simulation cell, molecular weight and number of chains 

The remaining parameters that define the amorphous eell are the number of ehains, 
built from to the monomer arrangements deseribed in seetion[2T](with eonfigura¬ 
tions determined by methods detailed in seetion lZ^ . the number of monomers in 
eaeh ehain (that determines moleeular weight), the size of the simulation eell, and 
the temperature. The simulation eell size ean be speeified direetly, as the lengths 
of a orthogonal eell along the three Cartesian direetions (in Angstroms) or as the 
density of a eubie box, whose dimensions are determined from the total mass and 
the density. 

2.3 Configuration building 

The prineipal polymer builder algorithm supported by PolymerModeler is the eon- 
tinuous eonfigurational biased direet Monte Carlo lUTll that effieiently samples tor¬ 
sions along the polymer baekbone and leads to struetures with reasonable statis- 
tieal weights. Chain eonformations are governed by eovalent interaetions assoei- 
ated with torsional angles along the baekbone and by non-bond interaetions that 
eause ehains to avoid overlaps and for affine groups to tend to eoalesee. These 
interaetions eause ehains to extend in good solvents and eoil in bad solvents. In 
addition to the eontinuous eonfigurational biased direet Monte Carlo algorithm 
PolymerModeler allows to build freely rotating ehain polymers, rod-like poly¬ 
mers and ehains built using self-avoiding walks for edueational and verifieation 
purposes. 
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2.3.1 Continuous Configurational Biased Direct Monte Carlo 

Due to their high stiffness, bond lengths and angles are fixed at their equilib¬ 
rium values throughout the building proeedure bond distanees; this is eommonly 
done in polymer builders. Thus, dihedral angles between baekbone atoms (0,) 
are the only degrees of freedom eonsidered during the initial strueture ereation. 
The first monomer is plaeed into the simulation eell at a random loeation (with 
uniform probability density) and with random orientation (with uniform angular 
distribution). Monomers are added sequentially with their torsion angles ehosen 
aeeording to the Boltzmann weight, i.e. the probability density of torsion i taking 
the angle 0* is proportional to exp [—E (0*) /kT] where E (0*) is the energy of 
the of the eorresponding eonfiguration, k is Boltzmann’s eonstant and T is the 
absolute temperature. The energy assoeiated with the eonfiguration is obtained by 
adding van der Waals interaetions and the eovalent torsional energy as deseribed 
in the following sub-seetions. 

The PolymerModeler builder uses a Metropolis (SSI Monte Carlo algorithm to 
explore and seleet torsion angle eonfigurations. Initial angles are assigned ran¬ 
domly for all flexible torsions in the new monomer and a Metropolis Monte Carlo 
Markov ehain is performed to determine the aetual torsion angle values aeeording 
to their Boltzmann weight. The Monte Carlo steps are as follows: 

• A new trial eonfiguration (torsion angles for all flexible dihedrals of the 
eurrent monomer) is ehosen by an unbiased, random modification of the 
current state, 

• The change in energy, AE, associated with the trial move is computed from 
the van der Waals interactions and covalent torsion terms (if specified); 

• The new configuration is accepted with probability one if /S.E < 0 (i.e. the 

Ag 

trial move decreases the energy of the system) and with probability e 
otherwise. 

This process is repeated for a user-specified number of steps, and the final angles 
are used to place the monomer. 

Covalent energy of backbone torsion angles The PolymerModeler amorphous 
builder allows users to specify the covalent torsional potential in two ways. The 
torsional energy as a function of angle can be provided numerically. Sample tor¬ 
sion potentials for sp^ and sp^ bonds according to the Dreiding force field are 
included in the tool. Alternatively, users can specify a probability as a function 
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of torsion angle; these probabilities are then converted into relative energies using 
Boltzmann weights. Energies for arbitrary atoms are obtained by interpolating 
between the tabular data. 

Non-bond interactions. Non-bond interactions in the continuous configura¬ 
tional biased direct Monte Carlo play a key role avoiding bad contacts between 
atoms. PolymerModeler computes the non-bond energy as the sum of Lennard- 
Jones pair interactions between each atom in the tail monomer and all surround¬ 
ing atoms, using Dreiding parameters. As is customary in molecular force fields, 
atoms in the same chain separated by fewer than four bonds are excluded from 
these non-bond interactions. The cutoff number of bonds for non-bond interac¬ 
tions is a configurable parameter. 

Computation of pair-wise interactions. A brute-force calculation of inter¬ 
atomic distance between the atoms of every new monomer with all previously built 
requires 1/2A^(A^ — 1) computations (N being the number of atoms) and becomes 
prohibitively expensive for relatively large systems. Thus, PolymerModeler uses 
the well-known domain decomposition approach to transform this computation 
into an 0{N) calculation. The simulation cell volume is divided via a 3D grid 
with spacing AL, close to but larger than the cutoff distance for van der Waals 
interactions for the Monte Carlo algorithm. As new monomers are added into the 
system atoms are placed into their corresponding cell; this involves 0{N) cal¬ 
culations. In order to compute inter-atomic distances to compute van der Waals 
energies or impose volume exclusions for an atom belonging to cell i only atoms 
in cell i and in its 26 nearest neighbors need to be checked. This calculations is 
independent of the simulation cell size for a given density and, thus, the algorithm 
is 0{N). 

2.4 Other Chain Builder Methods 

In addition to the configurational biased direct Monte Carlo, PolymerModeler also 
allows users to build polymer according to the methods described in the following 
paragraphs which are included for educations and verification purposes alone. 

Freely rotating chains. In this case all torsional angles have equal probabilities 
regardless of the covalent or van der Waals energies involved. 

Rod-like polymer chains can be created by specifying a single torsional angle 
to be imposed on all torsional degrees of freedom. 

Self-avoiding chains with a hard cutoff can also be created using excluded vol¬ 
umes. Excluded volumes are enforced by rejecting any configuration that places 
any 2 atoms within a cutoff distance. The attrition problems introduced by this 
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approach are well known [[Ml [35l [36l [37]| and we support it for edueational and 
verifieation purposes. 

2.5 Amorphous builder results 

This seetion demonstrates various eapabilities of the PolymerModeler amorphous 
builder ineluding output funetions and verifies its implementation by eomparing 
simulation results with analytieal solutions for known eases and well-established 
solutions. 

2.5.1 Freely rotating chains 

Figure shows the mean end-to-end length as a funetion of number of baek- 
bone bonds n obtained for 200 freely rotating ehains of isotaetie PMMA with 500 
monomers in eaeh ehain. Eaeh ehain was built independently of the others, with 
only one ehain in the simulation eell at any given time. The PolymerModeler re¬ 
sult is eompared with the exaet solution that in this ease is known analytieally as 
it eorresponds to a simple random walk [l24ll : 


{R) = l 


1 -I- cos9 


n 


1 — cosd 


( 1 ) 


with a C-C baekbone bond length I = 1.53A and baekbone bond angle 9 = 68°. 
The simulation result agrees with the analytieal solution. Figure [3] shows the dis¬ 
tribution of torsion angle values for the same simulation. Torsions ean take any 
value with equal probability in a freely rotating system and the eode produees the 
expeeted result. 

The distribution of end-to-end ehain lengths for the same PMMA system is 
shown in Figure HI The expeeted distribution Il2^ . the probability of finding the 
end of a ehain with length {R) in a spherieal shell with volume AirR'^dR, is also 
shown: 


P(R, {R))4nR‘dR = 4w ® 

with (R) ealeulated using ([1]). The simulation and theoretieal result agree well, 
and using a larger sample of ehains would improve the agreement. 
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Figure 2: End-to-end length of 200 ehains of freely rotating isotaetie PMMA with 
500 monomers eaeh 


2.5.2 Chains obtained via Continuous Configurational Biased Direct Monte 
Carlo 

In this seetion we present the results of building the same system as in the pre¬ 
vious seetion, 200 ehains of isotaetie PMMA with 500 monomers in eaeh ehain, 
but using the eontinuous eonfigurational bias Monte Carlo approaeh eonsidering 
Lennard-Jones interaetions between the atoms in the ehains (no eovalent interae- 
tions in this example). The eutoff range for the interaetions is 4 A. Eaeh ehain is 
built independently, as in the freely rotating ease, to faeilitate eomparisons with 
analytieal results. The Metropolis Monte Carlo sampling deseribed in seetion [231 
seleet torsions eonsidering van der Waals interaetions whieh avoids bad eontaets, 
produeing an effeetively self-avoiding ehain with a soft eore. 

Eigure [5] shows the end-to-end ehain length of the ehains in this system as a 
funetion of the number of baekbone bonds, n. The theoretieal trend lf24l in end-to- 
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Figure 3: Torsion angle histogram for 200 ehains of freely rotating isotaetie 
PMMA with 500 monomers eaeh 


end length as a funetion of n is similar to the random walk ease Q but inereases 

faster than the square root of the number of monomers. Self avoid walks show a 

^ 0.6 

dependeney that is also shown in the figure, demonstrating good agreement 
between PolymerModeler and theory. As expeeted, the short-range repulsive in- 
teraetions eause the ehains to swell with respeet to the freely rotating example, 
leading to longer ehains. 

The interaetions also ehange the distribution of torsion angles, shown in Figure 
Torsion angles between 90 and 270 degrees are more energetieally favorable. 
The distribution of end-to-end ehain lengths, shown in Figure |7l follows the 
trend of the expeeted distribution [(24| for self-avoiding ehains, 

P(x)inx^dx = p, 

whieh is the probability of finding the end of a self-avoiding ehain of length x = 
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Figure 4: End-to-end length histogram for 200 chains of freely rotating isotactic 
PMMA with 500 monomers each 


R/VW) in volume Anx'^dx. 


3 Structure relaxation and molecular dynamics 

After the desired atomistic structure is built using the configurational biased di¬ 
rect Monte Carlo method it needs to be relaxed to remove close contacts. The tool 
performs a series of relaxations involving scaled-down van der Waals interactions 
that accomplish this task and also can ameliorate entanglement problems. Fol¬ 
lowing these relaxations users can perform MD simulations to predict a variety 
of polymer properties including transport, glass transition temperature, and me¬ 
chanical response. PolymerModeler allows users to relax the structures generated 
and perform MD simulations with the Dreiding potential and the LAMMPS sim¬ 
ulation code. After the simulation is performed all the key data is displayed back 
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Figure 5: End-to-end length of 200 chains of self-avoiding isotactic PMMA with 
500 monomers each 


in the user’s web browser for analysis or download. This section describes the 
generation of an energy expression to compute total energies, stress, and forces 
that enable relaxation and dynamics followed by the methods to relax structures 
and options for MD simulations on the resulting atomistic models. 

3.1 Interatomic potential 

A general expression of the total potential energy in molecular force fields is ob¬ 
tained as a sum of energies due to valence or bonded interactions and non-bonded 
interactions: 

where Ur represents bond stretch interactions, Uq bond angle bending, rep¬ 
resents dihedral angle contributions, improper (out-of-plane) torsions, U^dw 
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Figure 6: Torsion angle histogram for 200 ehains of self-avoiding isotaetie PMMA 
with 500 monomers eaeh 


eorresponds to nonbonded van der Waals interaetions, and f4 is eleetrostatie in- 
teraetions originating from partial atomie eharges. A variety of foree fields are 
available for polymer simulations with speeifie parameterizations and funetional 
forms. PolymerModeler uses the DREIDING [|^ foree field, due to its wide 
applieability despite the relatively small number of parameters and beeause it is 
well validated for polymer ehemistry. Partial atomie ehargers are obtained with 
the Gasteiger approaeh |[^|4^ . The following subseetions deseribe the Polymer- 
Modeler implementation of the energy expression to moleeular meehanies and 
dynamies simulations. 

3.2 Atom types 

Dreiding[[38l and other moleeular foree fields use atom types to distinguish be¬ 
tween like atoms in different bonding environments and setup the energy press 
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Figure 7: End-to-end length histogram for 200 ehains of self-avoiding isotaetie 
PMMA with 500 monomers eaeh 


ion in Eq. HI For example, the bond streteh, angle bending and torsional po¬ 
tentials are different for sp^ and sp^ C atoms. Dreiding uses four eharaeters to 
speeify atom types, the first two being the element name (adding for elements 
represented by a single eharaeter). The third eharaeter denotes hybridization or 
loeal bonding (2 means sp^, 3 means sp^, and R denotes resonanee) and the fourth 
is used for atoms with implieit hydrogens. The atom type is determined by the 
element and the eonneetivity of a given atom and PolymerModeler automatieally 
types atoms in the monomers inputs by the users. If desired, users ean speeify 
Dreiding atom types when speeifying monomers using the z-matrix format. 

3.2.1 Covalent interactions 

Covalent energies in Dreiding inelude bond streteh, angle bending, torsions in¬ 
volving dihedral angles, and improper torsions. Identifying the terms eorrespond- 
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ing to these interactions require bond connectivities. PolymerModeler obtains 
connectivity information from the monomer z-matrix and checks for additional 
bonds with a distance-based criterion using atomic covalent radii. 

The procedure for generating covalent interaction is briefly described as fol¬ 
lows: 


• The list of bond stretch terms is obtained from the builder and used to gen¬ 
erate the energy expression terms as described below. 

• Bond angles are identified by listing bond pairs with a common atom central 
atom; 

• Dihedral angles ijkl are obtained for each bond jk by identifying atoms 
bonded to each atom in the central bond, based on bond information and 
atom linkage. 

• Improper torsions are identified based on bond information and atom link¬ 
age. An improper term ijkl is identified when three atoms jkl are bonded 
to a central one i. 

Once all individual terms are identified PolymerModeler assigns Dreiding 
force field parameters to each. Harmonic bonds are used with parameters de¬ 
pending on the two atoms making a bond. Cosine-harmonic angle terms are used 
with parameters depending only on the type of central atom. The parameters de¬ 
scribing torsion interactions in Dreiding depend on the hybridization of the central 
two atoms and the nature of these atoms is used to setup the energy expressions. 
Exceptions to this rule involve some torsions where the third atom is not an sp^ 
center and in those cases the force field types of three atoms are used. The pa¬ 
rameters for improper torsions are straightforward since they depend solely on the 
central atom of the improper term. 

3.2.2 Electrostatic interactions 

Partial atomic charges in the PolymerModeler tool are obtained for all atoms 
via the iterative algorithm of Gasteiger and Marsili |I3S |40l| that uses a partial 
equalization of orbital electronegativity. These charges are calculated in a self- 
consistent manner from the ionization potentials and electron affinities of the neu¬ 
tral atoms and of their corresponding single charge cations, and they depend only 
on the connectivity of the atoms involved. A total charge of zero is maintained 
for each chain. A refinement of the charge calculation was introduced in Ref. 
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|[4n for structures with resonant vr bonds, sueh as aromatie rings. This is not yet 
implemented in PolymerModeler. 

Users ean ehoose to deseribe long-range eleetrostaties via the partiele-partiele 
partiele-mesh (pppm), the reeommended method, or the Ewald method as im¬ 
plemented in LAMMPS. Alternative, a real-spaee eutoff ean be used to deseribe 
eleetrostaties. 

3.3 van der Waals interactions 

Both Lennard-Jones and exponential-6 (Buekingham) funetions ean be used for 
van der Waals interaetions in Dreding, and in PolymerModeler MD simulations 
are performed using the exponential-6 funetions sinee they provide a more aeeu- 
rate deseription |[38l 1211 l42ll . Lennard-Jones van der Waals funetions are used, 
however, for the initial relaxation of the strueture ereated by the builder sinee, 
unlike the exponential-6 funetion, it leads to repulsive interactions even for short 
interatomie distanees. 

3.4 Structural relaxation 

A system ereated by the amorphous builder must be relaxed before MD ealeula- 
tions ean be performed to remove close eontacts that result in very large forees. 
Generalizing the method of Theodorou and Suter [|43l, the PolymerModeler tool 
relaxes a system in multiple, sueeessive eonjugate gradient minimizations of the 
total energy. Sealed Dreiding Lennard-Jones parameters are used for this relax¬ 
ation. Eaeh minimization seales both the energy and distance parameters. The 
use of sealed van der Waals interaetions is motivated by the reduction in entangle¬ 
ments as ehains ean slide past one another as is done in state-of-the-art methods 
to relax amorphous polymers [[T9l . In addition, a gradual applieation of non-bond 
interaetions has been shown to avoid unphysieal ehanges in ehain eorrelation fune- 
tions Il44l as will be diseussed in more detail in Seetion [T5l 

Eigure[8]shows the energy of a first stage minimization, using LAMMPS, of an 
ataetie PMMA ehain of 96 monomers, built with the PolymerModeler amorphous 
builder at a density of 0.7 g/em^. Energies eorresponding to steps 4 through 500 of 
the minimization are shown; the first 3 steps are omitted for seale, as the energies 
of those steps are signifieantly larger. The relaxed system is then ready for MD, 
using unsealed Dreiding exponential-6 van der Waals interaetions. 
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Figure 8: Energy minimzation in LAMMPS of a single atactic PMMA chain with 
96 monomers 


3.5 Molecular dynamics 

After building and relaxing the structure is ready for MD simulations and prop¬ 
erty calculations. PolymerModeler provides users a large degree of control over 
the MD simulations that can be automatically performed within the tool using 
nanoHUB.org computational resources; users can also use PolymerModeler to 
generate LAMMPS input files and perform the simulations elsewhere. 

The PolymerModeler tool presents the user with a choice of MD ensembles, 
including NVT, NPT, and NVE. The following options are also configurable: 

• Time step 

• Number of steps 

• Temperature 
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• Temperature change per time step 

• System pressure 

The user may also choose to deform the system by specifying the engineering 
strain rate, true strain rate, or applied stress along each Cartesian axis. This is 
useful to characterize the mechanical response of the structures created via non¬ 
equilibrium MD simulations. 

Other MD options include the frequency of outputs, plotting the density as a 
function of time, and an initial thermalization prior to any MD steps. Furthermore, 
PolymerModeler presents the option to run LAMMPS in serial mode on a single 
processor core or in parallel up to 64 cores on nanoHUB HPC resources. The tool 
displays the progress of the MD simulation in the web browser window. 

After the simulation is performed the results are displayed graphically in the 
user’s web browser. These fully interactive outputs include: 

• Structure produced by the amorphous builder, the input to the MD simula¬ 
tion 

• An animation of two monomers rotating about the backbone bond between 
them 

• LAMMPS input file (command script) 

• LAMMPS data file 

• Plot of system energy vs. conjugate gradient step during minimization prior 
to MD 

• Plot of total and potential energies vs. time during MD simulation 

• Plot of kinetic energy vs. time during MD simulation 

• Plot of temperature vs. time during MD simulation 

• Diagonal components of the stress tensor vs. time during MD simulation 

• Off-diagonal (shear) components of the stress tensor vs. time during MD 
simulation 

• Animation of atoms during relaxation and MD simulation 
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Output log containing all text output from the amorphous builder and LAMMPS 


All outputs ean be downloaded as an image or as eomma-separated text, in the 
ease of plots. 



Figure 9: Time evolution of density for 5 PMMA ehains (1 ataetie, 4 syndiotaetie) 


Figure |9] shows a eomparison of two systems, eaeh with 5 PMMA ehains with 
96 monomers in eaeh ehain; 1 ehain was ataetie and 4 were syndiotaetie. One sys¬ 
tem was built at an initially high density with the eommereial paekage MAPS, and 
the other was built at an initially low density with the PolymerModeler amorphous 
builder. The built systems were then relaxed using LAMMPS (separate from the 
PolymerModeler tool) in NPT simulations at 600K for 500 ps with a time step of 
1 fs. The trend in density over time is very similar in both eases, indieating that 
the system built by PolymerModeler is as realistie as the system built by MAPS. 
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Using PolymerModeler to predict the structure and 
properties of amorphous polymers 

4.1 Molecular structure of polyethylene 

The molecular relaxation timescales of polymer melts (which increase with molec¬ 
ular weight as for high molecular weight cases HSl l are prohibitively long for 

molecular dynamics simulations in most cases of interest. Thus the initial molecu¬ 
lar structure and its relaxation is critical for accurate MD simulations. [|4^l44ll The 
mean square displacement (MSD) between monomers as a function of their index 
in a polymer chain was found to be a useful descriptor of the polymer structure. 
Auhl et al. Il44l showed that a simple relaxation of as-built polymer samples via 
energy minimization leads to unphysical molecular structure changes that relax 
over very long timescales and, consequently, affect predictions. These changes 
are reflected as significant deviations in the MSD for short and intermediate dis¬ 
tances. They found that a slow push-ojf approach, where intermolecular forces 
are turned on gradually, minimizes these problems and leads to well-equilibrated 
molecular structures. 

PolymerModeler allows multi-step relaxations where non-bond interactions 
are turned on gradually and we subjected our model to the same tests as in Ref. 
[|44l . Figure[l0]shows the mean square internal distances (MSD) between monomers 
for chains of polyethylene built using various algorithms, for comparison with 
Ref. Il44l . As expected, building the chains using a naive, freely-rotating model 
without avoiding close contacts leads to the smallest mean square displacements; 
adding Lennard Jones interactions during build increases the MSD and using the 
configuration biased Monte Carlo approach including van der Waals interactions 
and the torsional barriers caused by covalent terms leads to an even higher MSD 
and a more accurate structure. For this last cases, the system was built with sp^ 
torsional potential from Dreiding (Equation 13 in reference [|38l) and Lennard- 
Jones interactions between all atoms to a cutoff of 6A. The as-built structure was 
then relaxed in a gradual, 10-step, procedure where van der Waals interactions are 
turned on slowly. At each step we performed 500 conjugate gradient steps and 
500 molecular dynamics steps at a temperature T=600 K. The properties of the 
resulting structure are shown as dashed lines in Fig. [TOl Finally, we performed an 
additional thermalization with MD for 50,000 steps under isothermal, isochoric 
conditions at 600 K. The MSD of the structure after the MD simulations is shown 
as a solid line in Fig. [TOl Consistent with the findings by Auhl et al. we find that 
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a gradual relaxation procedure results in structures that maintain their statistical 
properties during relaxation and dynamics. Such structures produced by the Poly- 
merModeler builder with sp^ torsions and Lennard-Jones interactions followed by 
a gradual introduction of non-bond interactions are expected to be good starting 
points for MD simulations. 

In addition to the MSD, the radial distribution function g(r) for chains con¬ 
structed by PolymerModeler demonstrates good equilibration and good agreement 
with neutron diffraction experiments. FiguredUshows the radial distribution func¬ 
tion, g(r), for backbone carbon atoms in 100 isotactic PMMA chains built at 0.6 
g/cm^, fully relaxed with the 10 steps of minimization, after 400 ps of MD at 
600K. The peaks corresponding to nearest neighbor distances are labeled. In¬ 
dicators at 4. 8A and 8.6A are included in figure [H] for comparison with Figure 
7(a) in Ref. [|47l . which shows a similar g(r) calculation for the backbone atoms of 
syndiotactic PMMA. The indicators correspond to intrachain peaks (4.8A) and in¬ 
terchain peaks (8.6A) identified in PMMA by neutron scattering data. Structures 
built and equilibrated by PolymerModeler compare well to these results. 

We encourage PolymerModeler users to perform similar tests of their struc¬ 
tures before production runs. We also note that the structures generated can used 
as starting points for advanced parallel-replica type approaches specifically de¬ 
signed to equilibrate long-chain amorphous structures [|T9l . 

4.2 Prediction of the glass transition temperature of polymthyl- 
methacrylate 

In this section we illustrate the prediction of the glass transition temperature (Tg) 
of polymthylmethacrylate (PMMA) using MD over a structure built with Polymer- 
ModelerWe created a molecular structure using the continuous configurational bi¬ 
ased direct Monte Carlo algorithm using torsional potentials corresponding to sp^ 
bonds and a cutoff of 6 A to compute van der Waals interactions. The simulation 
cell contains 40 chains, each 90-monomers long, for a total of 54,080 atoms. The 
structure is then relaxed using a 3 phase minimization and thermalized for 100 ps 
under NPT conditions at 600 K. In order to compute Tg the relaxed structure is 
cooled down to room temperature under constant pressure conditions. We control 
the temperature and pressure of the system using Nose-Hoover thermostat and 
barostat with coupling timescales of 0.1 ps and 1 ps respectively. The system is 
cooled down in continuous increments at a rate of 0.1 K/ps. These simulations 
take approximately 70 hours when run using 24 cores. 
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Monomer separation, n 


Figure 10: Mean square internal distanees between monomers for 100 ehains of 
polyethylene built at 1.0 g/em^. Built systems were then relaxed and thermalized. 


Figurefl^shows the resulting density as a funetion of temperature during eool- 
ing. The thermal eontraetion that results from the anharmonieity in the interaetion 
potentials is elearly seen. More importantly for our study is the ehange is slope 
that ean be observed around 420 K. This denotes the glass transition where the 
system transforms from a melt to a glassy solid. The molten state above Tg has a 
higher thermal expansion eoeffieient that the solid below it. The predieted value 
of Tg is in good agreement with experimental values 388 K [|48l and 385 K [l49ll . 
When eomparing MD and experimental glass transition temperatures the effeet 
of eooling/heating rate should be eonsidered. The glass transition temperature in- 
ereases with inereasing heating/eooling rate and this inerease is approximately 3K 
per order of magnitude inerease in rate llSOllSTlISll . 
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Figure 11: Radial distribution function for PMMA backbone carbon atoms for 
100 chains of isotactic PMMA. 


5 Conclusions 

The PolymerModeler tool at nanoHUB.org features an amorphous builder that 
allows customized specification of monomers, including multiple monomers and 
monomer patterns to support stereochemistry. Chain conformations are deter¬ 
mined by specifying the torsion angle probabilities of backbone atoms. The 
well-known cases of freely rotating and self-avoiding chains are accurately re¬ 
produced by the builder. Structural properties of systems built by the Polymer- 
Modeler builder compare well with systems built by commercial packages and 
current state-of-the-art techniques. 

After a system of polymer chains is built it may be relaxed, thermalized, 
deformed, and otherwise simulated using LAMMPS directly from the Polymer- 
Modeler tool. Large LAMMPS simulations run on HPC resources affiliated with 
nanoHUB. All simulations are launched directly from a standard web browser. 
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Figure 12: Glass transition while cooling 54,080 PMMA atoms. 


and all results are displayed in the browser. No additional software, installation, 
or remote accounts are necessary to use this free tool. The structures produced 
by PolymerModeler can also be used as starting points for additional relaxations 
using special-purpose algorithms IfT^ . 

Users can submit support tickets and ask questions in nanoHUB forums. The 
PolymerModeler tool, like all nanoHUB tools, features a wish list of requested 
features. Future development is guided in part by requests from users via the wish 
list. Some of the features we would like to incorporate in PolymerModeler are: 

• Partial charge calculations using electronegativity equalization methods. 

• Implementing force fields other than Dreiding for MD simulations. 

• Gasteiger charges for resonant vr bonds 

We believe the combination of a powerful and flexible amorphous builder, the 
ability to perform LAMMPS-based MD simulations, the convenience of running 
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simulations interactively through a web browser, and the user feedback options 
make PolymerModeler a very compelling, free alternative to commercial polymer 
simulation packages. Combining power and ease of use, we foresee PolymerMod¬ 
eler will be used both for research and education. 
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